Perturbation Trasmission and Graph visualization¶

In [ ]:
# Import necessary libraries
import networkx as nx  # For creating and manipulating graphs
import numpy as np  # For numerical operations
import matplotlib.pyplot as plt  # For plotting
from matplotlib.patches import Patch  # For creating custom legend patches

# Set the plotting style to 'ggplot'
plt.style.use('ggplot')

# Set the default figure size to [10, 10]
plt.rcParams['figure.figsize'] = [10, 10]

SIRS model simulation¶

In [ ]:
def dW(delta_t: float) -> float:
    """Sample a random number at each call."""
    return np.random.normal(loc=0.0, scale=np.sqrt(delta_t))

def Euler_Maruyama_method(parameters):
    """
    Perform Euler-Maruyama method to simulate an epidemiological model.

    Parameters:
    - parameters (dict): Dictionary containing model parameters.
        Expected keys:
        't_in': Start time of simulation.
        't_end': End time of simulation.
        'N': Number of time steps.
        'mu': Natural birth and death rate.
        'b0': Transmission parameter.
        'b1': Transmission parameter.
        'phi': Transmission parameter.
        'gamma': Recovery rate.
        'ni': Disease-induced death rate.
        'alpha': Scaling factor for stochastic term.
        'S_in': Initial susceptible population.
        'I_in': Initial infectious population.
        'R_in': Initial recovered population.

    Returns:
    - TS (numpy.ndarray): Array of time steps.
    - Ss (numpy.ndarray): Array of susceptible population over time.
    - Is (numpy.ndarray): Array of infectious population over time.
    - Rs (numpy.ndarray): Array of recovered population over time.
    """
    t_in = parameters['t_in']
    t_end = parameters['t_end']
    N = parameters['N']
    mu = parameters['mu']
    b0 = parameters['b0']
    b1 = parameters['b1']
    phi = parameters['phi']
    gamma = parameters['gamma']
    ni = parameters['ni']
    alpha = parameters['alpha']
    S_in = parameters['S_in']
    I_in = parameters['I_in']
    R_in = parameters['R_in']

    # Calculate time step size
    dt = float((t_end - t_in) / N)
    # Create array of time steps
    TS = np.arange(t_in, t_end + dt, dt)
    assert TS.size == N + 1

    # Initialize arrays for S, I, R populations
    Ss = np.zeros(TS.size)
    Is = np.zeros(TS.size)
    Rs = np.zeros(TS.size)

    # Set initial values
    Ss[0] = S_in
    Is[0] = I_in
    Rs[0] = R_in

    # Perform Euler-Maruyama method to simulate SIR model
    for i in range(1, TS.size):
        t = t_in + (i - 1) * dt
        S = Ss[i - 1]
        I = Is[i - 1]
        R = Rs[i - 1]

        # Calculate effective transmission rate with stochastic term
        b0_tilde = b0 + alpha * dW(dt)
        beta = b0_tilde * (1 + b1 * np.cos(2 * np.pi * t + phi))

        # Calculate stochastic term
        dW_dt = dW(dt)

        # Update S, I, R populations using differential equations
        Ss[i] = S + ((mu - mu * S - beta * S * I + gamma * R) * dt - (alpha / b0_tilde) * beta * S * I * dW_dt)
        Is[i] = I + ((beta * S * I - ni * I - mu * I) * dt - (alpha / b0_tilde) * beta * S * I * dW_dt)
        Rs[i] = R + (ni * I - mu * R - gamma * R) * dt

    return TS, Ss, Is, Rs

# Example parameters dictionary
parameters = {
    't_in': 0,
    't_end': 1,
    'N': 5000,
    'mu': 0.009,
    'b0': 36.4,
    'b1': 0.38,
    'phi': 1.07,
    'gamma': 1.8,
    'ni': 36,
    'alpha': 0.25,
    'S_in': 0.998,
    'I_in': 0.002,
    'R_in': 0.0
}

# Run simulation using Euler-Maruyama method
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)

# Calculate basic reproduction number R0
def calculate_R0(parameters):
    """
    Calculate the basic reproduction number R0 for an SIR-type epidemiological model.

    Parameters:
    - parameters (dict): Dictionary containing model parameters.
        Expected keys:
        'mu': Natural birth and death rate.
        'gamma': Recovery rate.
        'ni': Disease-induced death rate.
        'b0': Transmission parameter.
        'phi': Transmission parameter.
        'b1': Transmission parameter.

    Returns:
    - R0 (float): Basic reproduction number.
    """
    mu = parameters['mu']
    gamma = parameters['gamma']
    ni = parameters['ni']
    beta = parameters['b0'] * parameters['phi'] * parameters['b1']
    R0 = beta / (mu + ni + gamma)
    return R0

# Calculate R0
R0 = calculate_R0(parameters)

# Print basic reproduction number (R0)
print(f"Basic Reproduction Number (R0): {round(R0, 1)}")

# Interpretation of R0
if round(R0) < 1:
    print("The infection will likely die out over time.")
elif round(R0) == 1:
    print("The infection will remain stable in the population but not cause an outbreak.")
else:
    print("The infection will likely spread and cause an epidemic.")

# Additional statistics
peak_infection = max(Is)
final_size = Ss[-1] + Is[-1] + Rs[-1]
print(f"Peak Infection: {peak_infection}")
print(f"Final Size of the Epidemic: {final_size}")
Basic Reproduction Number (R0): 0.4
The infection will likely die out over time.
Peak Infection: 0.0030332320384369007
Final Size of the Epidemic: 1.0009432020818723

SIRS PLOTS¶

In [ ]:
# Function to create a plot
def plot_population(TS, data, label, color):
    plt.figure(figsize=(10, 6))
    plt.plot(TS, data, label=label, color=color)
    plt.xlabel('Time')
    plt.ylabel('Proportion')
    plt.title(f'{label} Population Over Time')
    plt.legend()
    plt.grid(True)
    plt.show()

# Plot Susceptible over time
plot_population(TS, Ss, 'Susceptible', 'blue')

# Plot Infected over time
plot_population(TS, Is, 'Infected', 'red')

# Plot Recovered over time
plot_population(TS, Rs, 'Recovered', 'green')

# Plot all three (S, I, R) on the same graph
plt.figure(figsize=(10, 6))
plt.plot(TS, Ss, label='Susceptible', color='blue')
plt.plot(TS, Is, label='Infected', color='red')
plt.plot(TS, Rs, label='Recovered', color='green')
plt.xlabel('Time')
plt.ylabel('Proportion')
plt.title('SIR Model Over Time')
plt.legend()
plt.grid(True)
plt.show()
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image

NetworkX¶

1. Small-World Networks (Watts-Strogatz Model)¶

Small-world networks have high clustering like regular networks but also have short average path lengths between nodes. The Watts-Strogatz model can be used to generate small-world networks.

In [ ]:
np.random.seed(42)

# Generate a Watts-Strogatz small-world network
n = 300
k = 4  # Each node is connected to k nearest neighbors
p = 0.15  # Probability of rewiring each edge
G = nx.watts_strogatz_graph(n, k, p)

# Ensure there is at least one initially infected node
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
while not any(state == 'I' for state in initial_states):
    initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
nx.set_node_attributes(G, dict(zip(G.nodes, initial_states)), 'state')

# Define colors for each state and create a custom legend
state_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
legend_elements = [Patch(facecolor=color, edgecolor=color, label=state) for state, color in state_colors.items()]

# Fixed layout for consistent node positions in plots
pos = nx.spring_layout(G)

# Function to update node states based on SIRS dynamics
def update_states(G, t, Is, Rs, Ss):
    for node in G.nodes:
        state = G.nodes[node]['state']
        if state == 'S' and np.random.random() < Is[t]:
            G.nodes[node]['state'] = 'I'
        elif state == 'I' and np.random.random() < Rs[t]:
            G.nodes[node]['state'] = 'R'
        elif state == 'R' and np.random.random() < Ss[t]:
            G.nodes[node]['state'] = 'S'

# Simulation and plotting
infection_alive = True  # Flag to track if infection is still spreading

# Euler-Maruyama method simulation
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)

# Determine plot intervals
plot_intervals = np.linspace(0, len(TS) - 1, 10).astype(int)

# Simulation and plotting
for t_index in plot_intervals:
    # Update node colors based on their states
    node_color = [state_colors[G.nodes[node]['state']] for node in G.nodes]

    # Draw the graph with updated node colors
    plt.figure(figsize=(10, 6))
    nx.draw(G, pos, with_labels=False, node_color=node_color, node_size=50)
    plt.legend(handles=legend_elements, loc='best')
    plt.title(f'Small-World Network (Watts-Strogatz Model) at t={TS[t_index]:.2f}')
    plt.show()

    # Determine the range for the next state update
    next_t_index = plot_intervals[np.searchsorted(plot_intervals, t_index) + 1] if t_index != plot_intervals[-1] else len(TS)
    
    # Update states of nodes in the specified range
    for t in range(t_index, next_t_index):
        update_states(G, t, Is, Rs, Ss)
        
        # Check if all nodes are in susceptible state ('S')
        if all(G.nodes[node]['state'] == 'S' for node in G.nodes):
            infection_alive = False

    # If infection has died out, print message and break out of loop
    if not infection_alive:
        print("The infection has died out.")
        break

# Calculate effective reproduction number Rt over time
Rt = calculate_R0(parameters) * Ss / len(G.nodes)
print(f"Effective Reproduction Number (Rt) at the beginning: {Rt[0]}")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
The infection has died out.
Effective Reproduction Number (Rt) at the beginning: 0.0013022154795595405

Various metrics to analyze the graph's properties¶

In [ ]:
# Degree distribution
degree_sequence = [d for n, d in G.degree()]

# Average degree
avg_degree = sum(degree_sequence) / len(degree_sequence)

print(f"Average degree: {avg_degree}")

# Clustering coefficient
avg_clustering = nx.average_clustering(G)

print(f"Average clustering coefficient: {avg_clustering}")

# Average shortest path length
avg_shortest_path_length = nx.average_shortest_path_length(G)

print(f"Average shortest path length: {avg_shortest_path_length}")

# Diameter (may need to handle exceptions if graph is not connected)
try:
    diameter = nx.diameter(G)
    print(f"Diameter: {diameter}")
except nx.NetworkXError:
    print("Graph is not connected.")

# Degree centrality
degree_centrality = nx.degree_centrality(G)
# Example: Print degree centrality of the first few nodes
print("\nDegree Centrality\n")
for node in list(G.nodes())[:5]:
    print(f"Node {node}: Degree Centrality = {degree_centrality[node]}")

# Betweenness centrality
betweenness_centrality = nx.betweenness_centrality(G)

# Example: Print betweenness centrality of the first few nodes
print("\nBetweenness Centrality\n")
for node in list(G.nodes())[:5]:
    print(f"Node {node}: Betweenness Centrality = {betweenness_centrality[node]}")

# Assortativity
assortativity = nx.degree_assortativity_coefficient(G)
print(f"Degree Assortativity: {assortativity}")


# Global efficiency
global_efficiency = nx.global_efficiency(G)
print(f"Global Efficiency: {global_efficiency}")


# Local efficiency
local_efficiency = nx.local_efficiency(G)
print(f"Local Efficiency: {local_efficiency}")

# Degree histogram
degree_sequence = [d for n, d in G.degree()]
plt.hist(degree_sequence, bins='auto', alpha=0.75)
plt.xlabel('Degree')
plt.ylabel('Frequency')
plt.title('Degree Distribution')
plt.show()

# Community detection using Louvain method
import community

partition = community.best_partition(G)
print(f"Number of communities detected: {len(set(partition.values()))}")

# Community detection using Louvain method
partition = community.best_partition(G)
print(f"Number of communities detected: {len(set(partition.values()))}")

# Modularity
modularity = community.modularity(partition, G)
print(f"Modularity: {modularity}")

# Visualization of communities (optional)
pos = nx.spring_layout(G)
plt.figure(figsize=(12, 12))
nx.draw_networkx_nodes(G, pos, node_size=20, cmap=plt.cm.RdYlBu, node_color=list(partition.values()))
nx.draw_networkx_edges(G, pos, alpha=0.4)
plt.title('Graph with Communities (Louvain method)')
plt.show()
Average degree: 4.0
Average clustering coefficient: 0.31297142857142857
Average shortest path length: 5.746570281124498
Diameter: 12

Degree Centrality

Node 0: Degree Centrality = 0.020080321285140562
Node 1: Degree Centrality = 0.01606425702811245
Node 2: Degree Centrality = 0.020080321285140562
Node 3: Degree Centrality = 0.028112449799196783
Node 4: Degree Centrality = 0.020080321285140562

Betweenness Centrality

Node 0: Betweenness Centrality = 0.02705176631312766
Node 1: Betweenness Centrality = 0.017062091061331035
Node 2: Betweenness Centrality = 0.028341754057983386
Node 3: Betweenness Centrality = 0.07175942440050946
Node 4: Betweenness Centrality = 0.01647276093760317
Degree Assortativity: 0.029907628017542705
Global Efficiency: 0.20521358507622173
Local Efficiency: 0.3944952380952381
No description has been provided for this image
Number of communities detected: 16
Number of communities detected: 15
Modularity: 0.745494
No description has been provided for this image

2. Random graph (Erdős-Rényi)¶

In [ ]:
np.random.seed(42)

# Create an Erdős-Rényi random graph
G = nx.erdos_renyi_graph(n=200, p=0.15)  # n: number of nodes, p: probability of edge creation

# Assign initial states to nodes based on S_in, I_in, R_in
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
for node, state in zip(G.nodes, initial_states):
    G.nodes[node]['state'] = state

# Simulate SIRS dynamics on the graph
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)

pos = nx.spring_layout(G)  # Fixed layout for all iterations

# Determine the indices of the time points to plot
plot_intervals = np.linspace(0, len(TS) - 1, 10).astype(int)

# Define colors for each state
node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}

# Create custom legend
legend_elements = [
    Patch(facecolor='blue', edgecolor='blue', label='Susceptible (S)'),
    Patch(facecolor='red', edgecolor='red', label='Infected (I)'),
    Patch(facecolor='green', edgecolor='green', label='Recovered (R)')
]

for t_index in plot_intervals:
    node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
    node_color = [node_colors[G.nodes[node]['state']] for node in G.nodes]

    # Draw the graph with updated node colors
    nx.draw(G, pos, with_labels=True, node_color=node_color)
    plt.legend(handles=legend_elements, loc='best')
    plt.title(f'Random graph (Erdős-Rényi) at t={TS[t_index]:.2f}')
    plt.show()

    

    # Update states of nodes based on SIRS simulation results up to the next plot interval
    for t in range(t_index, len(TS) if t_index == plot_intervals[-1] else plot_intervals[plot_intervals.tolist().index(t_index) + 1]):
        for node in G.nodes:
            if G.nodes[node]['state'] == 'S':
                if np.random.random() < Is[t]:  # Infection dynamics
                    G.nodes[node]['state'] = 'I'
            elif G.nodes[node]['state'] == 'I':
                if np.random.random() < Rs[t]:  # Recovery dynamics
                    G.nodes[node]['state'] = 'R'
            elif G.nodes[node]['state'] == 'R':
                if np.random.random() < Ss[t]:  # Loss of immunity dynamics
                    G.nodes[node]['state'] = 'S'

# Calculate effective reproduction number Rt over time
Rt = R0 * Ss / len(G.nodes)
print(f"Effective Reproduction Number (Rt) at the beginning: {Rt[0]}")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
Effective Reproduction Number (Rt) at the beginning: 0.005045444444444444

3. Scale-Free Networks (Barabási-Albert Model)¶

Scale-free networks are characterized by a small number of highly connected nodes (hubs) and many nodes with only a few connections. The Barabási-Albert model is a popular method to generate scale-free networks.

In [ ]:
np.random.seed(42)

# Generate a Barabási-Albert scale-free network
n = 200
m = 3  # Number of edges to attach from a new node to existing nodes
G = nx.barabasi_albert_graph(n, m)

# Assign initial states to nodes based on S_in, I_in, R_in
initial_states = np.random.choice(['S', 'I', 'R'], size=len(G.nodes), p=[parameters['S_in'], parameters['I_in'], parameters['R_in']])
for node, state in zip(G.nodes, initial_states):
    G.nodes[node]['state'] = state

# Simulate SIRS dynamics on the graph
TS, Ss, Is, Rs = Euler_Maruyama_method(parameters)

# Simulate SIRS dynamics on the graph
pos = nx.spring_layout(G)  # Fixed layout for all iterations

# Determine the indices of the time points to plot
plot_intervals = np.linspace(0, len(TS) - 1, 10).astype(int)

# Define colors for each state
node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}

# Create custom legend
legend_elements = [
    Patch(facecolor='blue', edgecolor='blue', label='Susceptible (S)'),
    Patch(facecolor='red', edgecolor='red', label='Infected (I)'),
    Patch(facecolor='green', edgecolor='green', label='Recovered (R)')
]

for t_index in plot_intervals:
    node_colors = {'S': 'blue', 'I': 'red', 'R': 'green'}
    node_color = [node_colors[G.nodes[node]['state']] for node in G.nodes]

    # Draw the graph with updated node colors
    nx.draw(G, pos, with_labels=True, node_color=node_color)
    plt.legend(handles=legend_elements, loc='best')
    plt.title(f'Scale-Free Networks (Barabási-Albert Model) at t={TS[t_index]:.2f}')
    plt.show()

    

    # Update states of nodes based on SIRS simulation results up to the next plot interval
    for t in range(t_index, len(TS) if t_index == plot_intervals[-1] else plot_intervals[plot_intervals.tolist().index(t_index) + 1]):
        for node in G.nodes:
            if G.nodes[node]['state'] == 'S':
                if np.random.random() < Is[t]:  # Infection dynamics
                    G.nodes[node]['state'] = 'I'
            elif G.nodes[node]['state'] == 'I':
                if np.random.random() < Rs[t]:  # Recovery dynamics
                    G.nodes[node]['state'] = 'R'
            elif G.nodes[node]['state'] == 'R':
                if np.random.random() < Ss[t]:  # Loss of immunity dynamics
                    G.nodes[node]['state'] = 'S'


# Calculate effective reproduction number Rt over time
Rt = R0 * Ss / len(G.nodes)
print(f"Effective Reproduction Number (Rt) at the beginning: {Rt[0]}")
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
Effective Reproduction Number (Rt) at the beginning: 0.005045444444444444